//
// File: computeperiodogram.cpp
//
// MATLAB Coder version            : 5.4
// C/C++ source code generated on  : 23-Apr-2025 16:26:18
//

// Include Files
#include "computeperiodogram.h"
#include "FFTImplementationCallback.h"
#include "freqDomainHRV_data.h"
#include "freqDomainHRV_rtwutil.h"
#include "psdfreqvec.h"
#include "rt_nonfinite.h"
#include "coder_array.h"
#include <cmath>
#include <string.h>

// Function Definitions
//
// Arguments    : const ::coder::array<double, 1U> &x
//                const ::coder::array<double, 1U> &win
//                double nfft
//                double Fs
//                ::coder::array<double, 1U> &Pxx
//                ::coder::array<double, 1U> &F
// Return Type  : void
//
namespace coder {
void computeperiodogram(const ::coder::array<double, 1U> &x,
                        const ::coder::array<double, 1U> &win, double nfft,
                        double Fs, ::coder::array<double, 1U> &Pxx,
                        ::coder::array<double, 1U> &F)
{
  array<creal_T, 1U> Xx;
  array<double, 2U> costab;
  array<double, 2U> sintab;
  array<double, 2U> sintabinv;
  array<double, 2U> wrappedData;
  array<double, 1U> b_xw;
  array<double, 1U> xw;
  double b_win;
  int acoef;
  int bcoef;
  int csz_idx_0;
  int loop_ub_tmp;
  if (std::isnan(Fs)) {
    Fs = 6.2831853071795862;
  }
  acoef = win.size(0);
  bcoef = x.size(0);
  if (acoef <= bcoef) {
    bcoef = acoef;
  }
  if (win.size(0) == 1) {
    csz_idx_0 = x.size(0);
  } else if (x.size(0) == 1) {
    csz_idx_0 = win.size(0);
  } else if (x.size(0) == win.size(0)) {
    csz_idx_0 = x.size(0);
  } else {
    csz_idx_0 = bcoef;
  }
  acoef = win.size(0);
  bcoef = x.size(0);
  if (acoef <= bcoef) {
    bcoef = acoef;
  }
  if (win.size(0) == 1) {
    bcoef = x.size(0);
  } else if (x.size(0) == 1) {
    bcoef = win.size(0);
  } else if (x.size(0) == win.size(0)) {
    bcoef = x.size(0);
  }
  xw.set_size(bcoef);
  if (csz_idx_0 != 0) {
    acoef = (x.size(0) != 1);
    bcoef = (win.size(0) != 1);
    csz_idx_0--;
    for (int k{0}; k <= csz_idx_0; k++) {
      xw[k] = x[acoef * k] * win[bcoef * k];
    }
  }
  loop_ub_tmp = static_cast<int>(nfft);
  b_xw.set_size(loop_ub_tmp);
  for (csz_idx_0 = 0; csz_idx_0 < loop_ub_tmp; csz_idx_0++) {
    b_xw[csz_idx_0] = 0.0;
  }
  if (xw.size(0) > nfft) {
    if (xw.size(0) == 1) {
      wrappedData.set_size(1, loop_ub_tmp);
      for (csz_idx_0 = 0; csz_idx_0 < loop_ub_tmp; csz_idx_0++) {
        wrappedData[csz_idx_0] = 0.0;
      }
    } else {
      wrappedData.set_size(loop_ub_tmp, 1);
      for (csz_idx_0 = 0; csz_idx_0 < loop_ub_tmp; csz_idx_0++) {
        wrappedData[csz_idx_0] = 0.0;
      }
    }
    bcoef = div_s32(xw.size(0), static_cast<int>(nfft));
    csz_idx_0 = bcoef * static_cast<int>(nfft);
    acoef = (xw.size(0) - csz_idx_0) - 1;
    for (int k{0}; k <= acoef; k++) {
      wrappedData[k] = xw[csz_idx_0 + k];
    }
    csz_idx_0 = acoef + 2;
    for (int k{csz_idx_0}; k <= loop_ub_tmp; k++) {
      wrappedData[k - 1] = 0.0;
    }
    for (acoef = 0; acoef < bcoef; acoef++) {
      csz_idx_0 = acoef * static_cast<int>(nfft);
      for (int k{0}; k < loop_ub_tmp; k++) {
        wrappedData[k] = wrappedData[k] + xw[csz_idx_0 + k];
      }
    }
    acoef = wrappedData.size(0) * wrappedData.size(1);
    for (csz_idx_0 = 0; csz_idx_0 < acoef; csz_idx_0++) {
      b_xw[csz_idx_0] = wrappedData[csz_idx_0];
    }
  } else {
    b_xw.set_size(xw.size(0));
    acoef = xw.size(0);
    for (csz_idx_0 = 0; csz_idx_0 < acoef; csz_idx_0++) {
      b_xw[csz_idx_0] = xw[csz_idx_0];
    }
  }
  if ((b_xw.size(0) == 0) || (static_cast<int>(nfft) == 0)) {
    Xx.set_size(loop_ub_tmp);
    for (csz_idx_0 = 0; csz_idx_0 < loop_ub_tmp; csz_idx_0++) {
      Xx[csz_idx_0].re = 0.0;
      Xx[csz_idx_0].im = 0.0;
    }
  } else {
    boolean_T useRadix2;
    useRadix2 =
        ((static_cast<int>(nfft) > 0) &&
         ((static_cast<int>(nfft) & (static_cast<int>(nfft) - 1)) == 0));
    internal::FFTImplementationCallback::get_algo_sizes(
        static_cast<int>(nfft), useRadix2, &acoef, &bcoef);
    internal::FFTImplementationCallback::generate_twiddle_tables(
        bcoef, useRadix2, costab, sintab, sintabinv);
    if (useRadix2) {
      Xx.set_size(loop_ub_tmp);
      if (static_cast<int>(nfft) > b_xw.size(0)) {
        Xx.set_size(loop_ub_tmp);
        for (csz_idx_0 = 0; csz_idx_0 < loop_ub_tmp; csz_idx_0++) {
          Xx[csz_idx_0].re = 0.0;
          Xx[csz_idx_0].im = 0.0;
        }
      }
      if (static_cast<int>(nfft) != 1) {
        internal::FFTImplementationCallback::doHalfLengthRadix2(
            b_xw, Xx, static_cast<int>(nfft), costab, sintab);
      } else {
        acoef = b_xw.size(0);
        if (acoef > 1) {
          acoef = 1;
        }
        Xx[0].re = b_xw[acoef - 1];
        Xx[0].im = 0.0;
      }
    } else {
      internal::FFTImplementationCallback::dobluesteinfft(
          b_xw, acoef, static_cast<int>(nfft), costab, sintab, sintabinv, Xx);
    }
  }
  psdfreqvec(nfft, Fs, F);
  b_win = 0.0;
  acoef = win.size(0);
  for (csz_idx_0 = 0; csz_idx_0 < acoef; csz_idx_0++) {
    b_win += win[csz_idx_0] * win[csz_idx_0];
  }
  Pxx.set_size(Xx.size(0));
  acoef = Xx.size(0);
  for (csz_idx_0 = 0; csz_idx_0 < acoef; csz_idx_0++) {
    double Xx_im_tmp;
    double Xx_re;
    double Xx_re_tmp;
    Xx_re_tmp = Xx[csz_idx_0].re;
    Xx_im_tmp = Xx[csz_idx_0].im;
    Xx_re = Xx_re_tmp * Xx_re_tmp - Xx_im_tmp * -Xx_im_tmp;
    if (Xx_re_tmp * -Xx_im_tmp + Xx_im_tmp * Xx_re_tmp == 0.0) {
      Xx_re /= b_win;
    } else if (Xx_re == 0.0) {
      Xx_re = 0.0;
    } else {
      Xx_re /= b_win;
    }
    Pxx[csz_idx_0] = Xx_re;
  }
}

} // namespace coder

//
// File trailer for computeperiodogram.cpp
//
// [EOF]
//
